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fvj around the final, equilibrium state. This approximation reproduces the expectation value 

C^^ of the boundary stress tensor with a 20% accuracy. Here we elaborate on these results 

iy-s and extend them to observables that are directly sensitive to the bulk interior, focusing 

_il for simplicity on the entropy production on the event horizon. We also consider next-to- 

^D leading-order corrections and show that the leading terms alone provide a better description 
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1 Introduction 

One of the most interesting recent developments in the gauge-gravity duality (holography) 
[2-4] is the study of thermalization processes in strongly coupled plasmas by means of their 
dual gravity description [1, 5-11]. 

Part of the motivation for these investigations comes from heavy ion collision (HIC) ex- 
periments at RHIC and LHC, in which the successful phenomenological description of data 
requires the applicability of viscous hydrodynamics with the quark-gluon plasma equation 
of state and a tiny shear viscosity r]/ s = 0{\/Att) about 1 fm/c or less after the collision 
[12]. While this thermalization time is puzzlingly small from the viewpoint of bottom-up 
weak-coupling thermalization, it lies within the ball-park suggested by holographic calcu- 
lations [1, 5-11]. Although the dynamics in a real HIC surely involves a combination of 
weakly and strongly coupled physics, the hope is that understanding the limit in which 
the physics is assumed to be strongly coupled at all length scales may help us bracket the 
real-world situation — see [13] for a review of applications of holography to the QGP. 



In the gauge-gravity duality the gravitational representation of the thermal state of a 
holographic quantum field theory in the planar limit and at strong ('t Hooft) coupling is 
a bulk black brane [14]. Thus thermalization on the gravity side is the process whereby a 
bulk geometry approaches the form of (a patch of) a static black brane and the other bulk 
fields relax to their form in the black brane background. Generically, this stage is preceded 
by one which one may dub 'hydrodynamization' whereby the system reaches a phase in 
which the system is still evolving but the dynamics is governed by hydrodynamic equations. 
On the gravity side this implies that the bulk geometry must have the form dictated by 
the so-called fluid-gravity duality [15] (see e.g. [16] for a review). As black brane metrics 
already depend on one coordinate and we want to allow for at least time dependence, the 
studies of holographic thermalization processes require solving Einstein's equations in two 
or more variables, which in most cases requires the use of numerical relativity techniques. 

Apart from toy-models based on the AdS-Vaidya geometry of infalling dust (see e.g. 
[17-20], but also [21]), the only existing approximation scheme to study holographic ther- 
malization processes is the amplitude expansion, in which one linearizes Einstein's equa- 
tions on top of the static black brane background. In this approximation the relaxation 
towards equilibrium is described by quasinormal modes with complex frequencies, whose 
imaginary parts lead to the damping of their amplitudes with time and hence to equilibra- 
tion. These modes were thought so far to be appropriate for the description of only the 
late-time approach to equilibrium, when deviations from equilibrium are sufficiently small 
in amplitude [22]. 

An indication that this assumption might be too restrictive comes from black hole 
mergers in asymptotically flat four-dimensional spacetime. There, in the so-called close- 
limit approximation, the Einstein's equations linearized on top of the flnal black hole predict 
rather accurately the pattern of gravitational radiation at infinity provided the initial data 
have a single horizon surrounding the merging black holes [23, 24]. This initial horizon, 
however, needs not to be a small perturbation of the final black hole for the close-limit 
approximation to work. 

These features, together with the observation that the AdS analogue of gravitational 
radiation at infinity is the expectation value of the boundary stress tensor, motivated 
recent work [1] in which a linear approximation was applied to a simple example of far-from- 
equilibrium gravitational dynamics in AdS spacetime. Specifically, this reference considered 
the isotropization of a large number of initially homogenous but anisotropic states in a 
Conformal Field Theory with a holographic description, or holographic Conformal Field 
Theory (hCFT) for short, with a large number of colors Nc and strong 't Hooft coupling 
[25] A = ^YM-^c- The outcome of this analysis was that the linearized Einstein's equations 
indeed reproduce well the dynamics of the one-point function of the stress tensor. For small 
perturbations this (expectedly) provided excellent approximations to the time evolution 
of the stress tensor. Surprisingly, however, even for initial states which were not a small 
perturbation of the final state, the linearized evolution still predicted the stress-tensor with 
an accuracy of about 20%. The purpose of the present paper is to explore the applicability 
of the linearized Einstein's equations further, in particular for observables that are directly 
sensitive to the interior bulk dynamics such as the entropy production. 



We start this article by reviewing the results of [1]. The key technical novelty of [1], 
compared to earlier works, was the ability to generate and follow the evolution of a large 
number (around 1000) of different non-equilibrium initial states. In the current article we 
provide a detailed explanation of how these initial states were obtained and evolved. We 
also comment on the relation between features of the initial states and the time it takes 
them to thermalize. 

Subsequently, we review the linearized gravity approach to holographic isotropization 
and comment on the relation between the form of initial data and the accuracy of the linear 
approximation. We also study in detail the decomposition of the solutions into quasinormal 
modes. According to the common lore, quasinormal modes can only describe the late-time 
approach to equilibrium, for which only the lowest-lying modes are relevant. As already 
anticipated in [1], our finding is that, provided a sufficient number of modes is considered, 
the description in terms of quasinormal modes can actually be extended all the way back 
to early times for which the system is very far from equilibrium. The correct description of 
the system at all times then results from the heavy interference of the quasinormal modes 
with one another. 

Finally, we discuss the corrections to the linearized Einstein's equations. For generic 
observables that are sensitive to the geometry deep in the bulk these corrections are crucial 
because they constitute the leading-order result. This is the case, for example, for the 
entropy production, on which we focus our attention. 

2 Holographic setup 

2.1 Holographic isotropization and the ansatz for the dual metric 

The holographic isotropization that we will study is the dynamics of a hCFT, in which 
a strongly coupled (3 + l)-dimensional plasma is homogeneous, but has a time-dependent 
pressure anisotropy leading to rich non-equilibrium physics [1, 5]. In this setup all the local 
operators apart from the stress tensor have vanishing expectation values, which constitutes 
the universal sector of any hCFT [26] . 

This example of a holographic thermalization process is, most likely, the simplest 
one, as due to homogeneity it does not involve hydrodynamic modes and thermalization 
is achieved directly via relaxation to a thermal state without being preceded by a long 
hydrodynamic phase. The latter observation is the primary reason for considering it here, 
and previously in [1], in the context of the close-limit approximation for anti-de Sitter 
spacetimes. 

As in [1], we will be interested in an unforced relaxation to the thermal state and 
we consider the situation in which there are no time-dependent sources pumping energy 
into our system, i.e. the boundary metric is flat. This is the crucial technical difference 
between our work and that of reference [5] , which considered the holographic isotropization 
process in which non-equilibrium states were obtained by placing the boundary theory in 
a time-dependent metric for a finite period of time. 

The form of the bulk metric is determined to a large extent by the ansatz for the 
boundary stress tensor. After imposing for simplicity rotational symmetry in two of the 



spacelike directions and taking into account the flatness of the boundary metric, the most 
general conserved and traceless stress tensor can be written as 



27r2' 



(V) = Td^diag £, V^{t), V^it), V^it) , (2.1) 



where N^ is the rank of the SU{N^) gauge group and £ is proportional to the energy 
density. The longitudinal Vi,{t) and transverse pressures VT{t) are then expressed in terms 
of a time-dependent anisotropy AP(t) through 

V^{t) = l£- lAV{t) , (2.2) 

VAt) = l£ + lAVit). (2.3) 

Note that the energy density in our setup is constant in time by virtue of the homogeneity 
assumption plus the conservation of the stress tensor. In this sense the energy density is 
part of the initial conditions. As the only possible static state with finite energy density is 
the isotropic and homogeneous plasma [27], the final state is known already from the start, 
without the need to solve any dynamical equation. This seems to be a rather non-generic 
feature of our setup, which we discuss at length in the conclusions section. 

The symmetries of the stress tensor - through the general near-boundary form of 
the metric in the Feffer man- Graham coordinates [28] - suggest the ansatz for the dual 
metric gab with gu, ^ll and g-j-T components being dynamical, i.e. obtained by solving the 
equations of motion.^ The freedom of defining what is meant by the bulk time t and the 
radial coordinate in AdS r has to be then (partly) fixed by imposing the form of the grr 
and gtr components. 

In our analysis it will be very convenient to follow [5] and express the dual metric gab 
in the generalized ingoing Eddington-Finkelstein coordinates 

ds^ = 2dtdr - Adt^ + T^e'^^dxl + T?e^d^l , (2.4) 

where A, T, and B are functions of time t and the radial coordinate r. With this ansatz 

grr = and gtr = 1. In the ingoing Eddington-Finkelstein coordinates, constant time slices 

are null hypersurfaces as radial ingoing null geodesies, by construction, penetrate the bulk 

instantaneously. 

The metric (2.4) has to solve the Einstein's equations with the negative cosmological 

constant 

1 6 

Rah - ^Rgab - J2 9ab = 0, (2.5) 

where in the following we set the radius of the AdS solution L to 1. The link between the 
form of the field theory stress tensor and the dual metric ansatz becomes clear after solving 



^It does not have to be necessarily the case and actually in the ADM formulation of general relativity 
it is the form of the lapse - the timelike warp factor (gtt) - and the shift (gta) that are fixed/imposed. For 
details in the context of AdS spacetimes see [9] . 



Einstein's equations in the near-boundary (large r) expansion 



2 a4 264 (t) 
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64(t) , dthit) 

B = ^:r + ^^ + • • • ' (2-6b) 

E = .-M^ + ..., (2.6c) 

and identifying nornializable modes with the components of the stress tensor using holo- 
graphic renormalization [28]. For the specific case of SU{N^) M = A super Yang-Mills 
theory this relation is 

£ = -2,ai/A and AP(t) = 364(i) • (2.7) 

The metric ansatz (2.4) leaves the residual diffeomorphisms freedom r — t- r + f{t). This 
freedom is fixed here by imposing the absence of the term proportional to r in the near- 
boundary expansion for the gu metric component (2.6a). 

In (2.6) we suppressed the near-boundary expansion at relatively low order, but it is 
important to stress that the expansion has infinitely many terms carrying arbitrarily high 
derivatives of the pressure anisotropy. This inevitably leads to a general conclusion that 
a state given by the form of the geometry on a constant time slice is (partly) specified by 
infinitely many derivatives of the dual stress tensor, in our case the pressure anisotropy. 

2.2 Thermalization criterium 

Although £ is constant in time, the physical temperature can only be assigned to the 
system once the equilibrium is reached. In this regime £ = 37r^T^/4 and the metric takes 
the form 



ds^ = 2dtdr -r^U- ^^^\ de + r'^dx^ , (2.^ 



or in terms of A, T, and B 



A = r^ [1- ^^^ ) , S = r and B = 0, (2.9) 



r^ 



and describes the AdS-Schwarzschild solution between the boundary and the future event 
horizon covering also the black brane interior. 

Although equilibration of a holographic system can be sampled with different field 
theory probes, including expectation values of local operators, two point functions, entan- 
glement entropy and Wilson loops, in our studies we will primarily focus on tracing the 
evolution of the one-point function of the stress tensor. There are two reasons for this. 
In the first place this is the quantity of interest if one wants to make a phenomenological 
contact with the fast applicability of hydrodynamics at RHIC and LHC. Secondly, after the 
stress tensor eventually settles down to its thermal value, the geometry becomes a patch 
of the AdS-Schwarzschild black brane and from this moment on there is no need to evolve 
the Einstein's equations further. 



We will hence define thermalization time as the isotropization time ti^o, i-e- the time 
after which the stress tensor anisotropy AV{t) remains small compared to the energy 
density and eventually decays to zero. In our calculations, as in [1], we adopt the following 
criterium for ti^^'- 

AV{t>t,,,)<0.1£, (2.10) 

but it is important to keep in mind that thermalization is never a clean-cut event and the 
threshold on the RHS of (2.10) can be always slightly raised or lowered without altering 
much the results. 

2.3 The event horizon and its entropy 

The event horizon is defined as the causal boundary of the black hole. It is a teleological 
object which can be located only after the black hole settles down to the state of permanent 
equilibrium. 

For the purposes of the current paper, we will be interested in the event horizon's area 
as an example of one easy-to-compute bulk observable that is directly sensitive to the form 
of the geometry in the deep IR. Although no precise definition of the entropy density exists 
in a truly far-from-equilibrium situation, the change in the area density of the event horizon 
provides a crude measure of the total entropy produced in the thermalization process. For 
this reason we will loosely refer to the area density of the event horizon as 'entropy density', 
but we emphasize from the start that this terminology is rigorously applicable only near 
equilibrium. We also clarify that we chose to focus on the event horizon, as opposed to 
e.g. the apparent horizon, because in many cases there was no apparent horizon on our 
initial-time slice. 

Due to the symmetries of our problem, the event horizon will be a hypersurface defined 

by 

r - r^nit) = 0, (2.11) 

with the normal vector being null 

rUt)-lA{t,rUt)) = 0. (2.12) 

The latter is the geodesic equation for the outgoing light ray and needs to be supplemented 
with the following condition to be imposed in the asymptotic future 

rEH(i) ^ vrT as t ^ oo . (2.13) 

In practical terms this condition implies that when the metric eventually approaches the 
form of the AdS-Schwarzschild black brane (2.8), Teh approaches the position of the event 
horizon of the static solution. 

The area of the event horizon gives rise to the following expression for the entropy 

SEH(t) = TT^iVc S {t, r^n{t)f , (2.14) 

zvr 

which is guaranteed to be a non-decreasing function of t. In (2.14) we implicitly as- 
sumed that the event horizon's area is mapped onto the boundary along ingoing null radial 
geodesies, i.e. along lines of constant t for the metric ansatz (2.4). This is a reasonable 
choice in our setup due to its high degree of symmetry. 



3 Far-from-equilibrium initial states and their nonlinear evolution 

3.1 Solving Einstein's equations as an initial- value problem 

As originally noted in [5], the Einstein's equations for the metric ansatz (2.4) can be neatly 
written as 

= S(S)' + 2S'i]-2S^ (3.1a) 

= ^{By + 1{T.'B + B't), (3.1b) 

= A" + 3-B'^ - 12S' S/S^ + 4 , (3.1c) 

= t + l{B^^-A't), (3.1d) 

= S" + iS'2s, (3.1e) 

where 

h' = drh and h = dth + \A drh (3.2) 

denote respectively derivatives along the ingoing and outgoing radial null geodesies. In the 
following we will be interested in solving the initial- value problem, i.e. given the geometry on 
the initial-time slice we want to obtain the evolution of the dual stress tensor by computing 
the bulk spacetime outside the event horizon. 

Not all equations among (3.1) are evolution equations, i.e. specify the form of the 
metric on a neighboring time slice. Equations (3. Id) and (3.1e) are constraints in the sense 
that the remaining components of the Einstein's equations can be shown to guarantee that 
they are obeyed provided (3. Id) holds at the boundary and (3.1e) holds on the initial-time 
slice [5]. 

The null character of the coordinate system (2.4) leads to a nested algorithm for solving 
the initial-value problem in which one uses as evolution equations (3.1a)-(3.1c) and at each 
time step one only needs to solve linear ordinary differential equations in r. The precise 
scheme that we will follow is a slight modification of the one originally introduced in [5] , 
and consists of the following steps: 

1 . we start with B as a function of r and the energy density £ (constant in our setup) ; 

2. the constraint equation (3.1e) allows us to solve for S as a function of r; 

3. we then solve (3.1a) for S, with £ being the integration constant; 

4. having B, S and S, we solve (3.1b) for B; 

5. with B, S, B and S at hand we can integrate (3.1c) for A; 

6. knowing B and A and using (3.2) we get dfB; 

7. we proceed to the next time step using a finite difference scheme (for details see 
section 3.3). 



In our setup the constraint (3. Id) is implemented in the near-boundary expansion, as it 
encodes the conservation of the stress tensor in the dual gauge theory. In order to monitor 
the accuracy of the numerical code we check the value of this constraint in the bulk when 
evaluated on the numerical solution (see also section 3.3). 

The algorithm outlined above needs to be supplemented with the initial conditions, 
the choice of which we discuss in the next subsection. 

3.2 Specifying initial states 

Gravity encodes dual initial states in the form of the geometry on a bulk initial-time slice. 
The conditions on the initial data arise from three sources: the constraint (3.1e), the near- 
boundary expansion (2.6) and bulk regularity. By the latter we mean that all possible 
singularities in the initial data must be hidden inside the event horizon. 

One way to obtain a non-equilibrium state while automatically satisfying the conditions 
above is to start with vacuum AdS and perturb it by turning on a non-normalizable mode 
of the bulk metric or some other bulk field for a finite period of time [5] . The alternative 
approach, that we adopt here and which was used also in [8, 9, 29], is to start with non- 
equilibrium states defined as solutions of the constraints on the initial-time slice without 
invoking the way in which a particular state was created. 

Equation (3.1e) imposes a constraint between the forms of B and S on the initial- 
time slice. Since B appears quadratically, we choose to specify the initial state through B 
and then use (3.1e) to solve for S. Note that this equation, together with the asymptotic 
behavior linear in r (2.6c), implies that S must necessarily be a convex function and hence 
that it must vanish for some r > 0, with the inequality being saturated only for vacuum 
AdS and the AdS-Schwarzschild black brane. As our coordinate frame is spanned by the 
ingoing radial null geodesies and S measures the transverse area of the congruence, S = 
implies reaching a caustic and hence the breakdown of our coordinate frame. 

For the successful evolution of the initial data specified by some B we thus need to 
make sure that the locus where S vanishes is hidden behind the event horizon on the 
initial-time slice. As the event horizon is a teleological object, this cannot be verified a 
priori - we need to try to run a simulation and when it is successful we know that the 
initial state we started with was legitimate. 

The contrary is not necessarily the case, as a caustic, a priori, is just a breakdown of a 
coordinate system. However, we verified numerically that in the neighborhood of a point 
where S vanishes we obtain very large curvatures. This suggests that this point must be 
hidden inside the event horizon.^ 

We thus see there is an interesting interplay between the choice of B and the choice of 
the (initial) energy density £. Both quantities, a priori, seem to be very much independent 
when it comes to specifying the initial state. If, however, the point where S vanishes 
corresponds to a genuine curvature singularity, which is the case for the AdS-Schwarzschild 
black brane and which our numerical studies also indicate, then there must be a minimal 

■^A possible caveat is that our numerics is not guaranteed to give trustable results inside the event 
horizon, as the physics might be influenced there by the artifacts of the inner cut-off of the grid. 



energy density E for which this singularity is still covered by the event horizon on the 
initial-time slice. An indication that this might be the case comes from noting that for 
the AdS-Schwarzschild black brane the radial position of the event horizon is proportional 
to (the fourth root of) the energy density E. Similarly, in our setup we can always put 
the final event horizon at a smaller value of r than the point at which S vanished on the 
initial-time slice by decreasing the energy density. This discussion suggests that the states 
for which the initial position of the event horizon is very close to the point where S = 
should be viewed as maximally far from equilibrium. 

In our setup, we have a freedom of preparing arbitrary initial conditions, i.e. we can 
specify i? as a function of r on the initial-time slice and E > 0, as long as B obeys the 
near-boundary expansion of the form (2.6b) and there are no naked singularities. We use 
this freedom to prepare and follow the evolution of states in which B has support very 
close to the boundary, very close to the horizon or spreads over a large range of the radial 
direction. In order to generate a large number of non-equilibrium initial states we followed 
the following procedure: 

1. without loss of generality we choose units such that a^ = —1, or equivalently i5 = |; 

2. we generate the initial i? as a ratio of two 10th order polynomials in r with random 
coefficients in the range (0, 1); 

3. we subtract from it a cubic expression so that the near-boundary expansion for B of 
the form (2.6b) is obeyed; 



4. the whole expression is then normalized so that the maximal value of the B between 

1. 

2' 



the boundary and the position of the final event horizon (r = 1) is -■ 



5. we then run a binary search algorithm to find the factor that B needs to be multiplied 
by such that the code is just stable, while storing successful runs. Typically, we repeat 
this step about 6 times per seed function generated in step 2. 

In this way we can generate states which are as far from equilibrium as our numerical 
code allows. In the end this means there is some sensitivity to the number of grid points, 
since increasing the number of grid points would improve the stability. 

Finally, it is interesting to note that a constraint of exactly the form (3.1e) also holds 
for metric ansatzes corresponding to a dual plasma expanding in one dimension [6, 7]. 
This implies that our discussion about the specification of the initial states, including the 
maximally far from equilibrium ones, also apply in these other setups. However, if we relax 
the assumption of a homogeneity in the transverse plane, then S is no longer forced to be 
convex and there might be bulk states which do not lead to caustics/apparent singularities 
in the way described above. 

3.3 Numerical implementation 

In the numerical implementation instead of the variable r we used its inverse 

z = -, (3.3) 

r 




Figure 1. The left plot shows the value of the normalized constraint £^{t — Q) as a function of the 
number of grid points n for the evolution of the initial profile B{z, ti^i) = l^^^- It is clear from 
the plot that our numerics converges exponentially with the number of grid points. The right plot 
shows the evolution of ^{t) as a function of time for n = 26 and one can see there that the constraint 
actually decreases with time. To achieve (,{t) < 10^^ one typically needs higher precision than the 
standard double precision computations offer. 

so that the AdS boundary is at z = 0. With the choice of units we made in our code, 
i.e. a4 = —1, the horizon of the final black brane is then located at z = 1. Note however 
that, for definiteness, all dimensionful quantities that we will provide will be specified in 
terms of the energy density or, equivalently, the temperature of the final black brane, which 
is the only scale at the moment of thermalization. 

For discretization we use the pseudo-spectral collocation method for the spatial grid 
in the z direction (see [30, 31] for accessible reviews of the spectral methods in the context 
of numerical GR). This allows us to maintain a very moderate grid, typically with 26 
points, without a significant loss of accuracy. Each of the evolution equations (3.1a)-(3.1c) 
is written in terms of redefined variables in which the three terms of lowest order in the 
near-boundary expansion are subtracted, and solved as a set of linear equations. 

The latter is a significant simplification, since it allows the use of standard procedures 
for solving linear equations. We implemented our numerical code in Mathematica and used 
the fourth order Adams-Bashforth method for evolution in time. The first five time steps 
were made with the use of low order Runge-Kutta methods. Typically, we used time steps 
of size 0.0025 ~ (# grid points)"^. 

As a way of monitoring the accuracy of our code, we used the normalized constraint (3. Id) 



^(t) = maXf, 



IB'^ 



lA't 



\t\ + lB'^^ + l\A't\ 



(3.4) 



fixed t 



The convergence of our code is then illustrated by Fig. 1, which shows typical plots of the 
maximum value of the normalized constraint ^{t). 

The last feature that needs to be discussed is the choice of the position of the inner 
boundary of the computational grid. Note that the simulation is well-defined only if the 
grid covers the entire portion of the spacetime outside the event horizon. Initially this is 
hard to predict, since the position of the event horizon depends on the future evolution. 
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Figure 2. The left plot shows B(z, t) for the initial profile (3.5), which is shown as a thick red 
line at i = 0. By construction, all time derivatives of the pressure anisotropy vanish at i = 0. The 
thick blue curve at z = shows the value of the gauge theory quantity A'P{t)/£. One sees clearly 
that the Taylor expansion around the initial time (predicting constant pressure anisotropy) does 
not converge when thermalization is achieved. 

The right plot, in which we start with B{t = 0, z) — (|f )^z^'', shows similar behavior. The initial 
disturbance, which is localized in the IR part of the geometry, propagates to the boundary in a 
time limited by causality. This creates the pressure anisotropy, which quickly relaxes back to zero. 

Therefore one typically focuses attention on the presence of the apparent horizon which, if 
out can be found, it is guaranteed to lie inside the black hole. However, quite frequently 
in our case there was no apparent horizon on the initial-time slice and therefore we used 
the following procedure. We first tried to run simulations with the radial cut-off put at 
z = 1.01, which is right below the late-time position of the event horizon. This often 
worked, and when it did not we reran the simulation with z = 1.07 as a cut-off. The 
latter point turned out to almost always lie past the event horizon. In this way we could 
successfully evolve a large number of initial states. 

3.4 Why the near-boundary expansion is not enough 

The near-boundary expansion of the metric contains information about all time derivatives 
of the pressure anisotropy at the initial time. This allows the construction of a Taylor-series 
expansion for the pressure anisotropy, which at least for some time should be very close to 
the exact expression. Previous studies in the context of boost-invariant flow suggest that 
the convergence radius of a series of this type is too small to see thermalization [9] . We 
can convince ourselves that this is the case also here by considering the initial profile of 
the form 



BU^) 



:£z' 



(3.5) 



According to the near-boundary expansion (2.6b) this profile at the initial time has a non- 
zero pressure anisotropy, but all its time derivatives vanish. The Taylor series expression 
predicts then a constant pressure anisotropy, which would lead to a singular spacetime [27]. 
Fig. 2 shows the numerical evolution of the profile (3.5) and proves that the radius 
of convergence of the Taylor expansion is indeed too small to see thermalization defined 
via (2.10). This result implies that the precise form of the metric Qab and the pressure 
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Figure 3. The left figure shows i? as a function of time and radial coordinate for the initial profile 
(3.6), which is shown as a thick red curve at t = and which is initially localized near z = ^Zh- 
The blue curve at the boundary (z = 0) depicts the pressure anisotropy as a function of time in the 
gauge theory. The right figure shows the Kretschmann scalar (with the value for an equilibrium 
black brane with the same energy density subtracted) as a function of time and radial coordinate 
for the same initial profile. One clearly sees on this plot the wave bouncing off the boundary and 
falling into the black brane. In the adopted generalized ingoing Eddington-Finkelstein coordinates 
this happens instantaneously. 

anisotropy AV{t) at transient times need to be obtained via an explicit solution of the 
initial-value problem on the gravity side. 

3.5 Evolution of a sample profile and expectations for thermalization times 

To get intuition about how the dynamics proceeds on the gravity side and to get acquainted 
with the features following from the choice of a foliation by null constant-time slices, it will 
be instructive to discuss in detail the dynamics of the following initial state 



B{t = 0,z) 
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tz exp 
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(3.6) 



where z^ = -^j^S^''^. As B is supported at intermediate values of z, naive intuition from 
the physics of linear wave equations would suggest that the wave packet splits into two: one 
propagating inwards and the other propagating outwards. The one propagating outwards 
is expected to eventually reach the boundary, bounce back and fall into the bulk. Both 
wave packets will be eventually absorbed by the event horizon (which is guaranteed to be 
present given that £ ^ 0) leading to the increase in its area. 

These expectations are confirmed by the outcome of the numerical simulation, as il- 
lustrated by Fig. 3, which depicts the bulk anisotropy (left plot) and the square of the 
Riemann tensor, the Kretschmann scalar (right plot). We can clearly see the rise in the 
curvature due to the outgoing wave packet as it approaches the boundary of AdS. Closer 
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inspection reveals also the presence of a wave packet resulting from the bouncing off the 
boundary of the outgoing packet. This wave packet, due to the null nature of our coordi- 
nate frame, propagates towards from the boundary to the horizon along lines of constant 
Eddington-Finkelstein time. Note also that this signal falls through the black brane event 
horizon without significant scattering. This feature persisted for other choices of initial 
states and seems to be related to the high degree of symmetry of our problem. 

It is interesting to note that the initial ingoing part of the wave packet seems to be 
mostly taken care of by the solution of the constraints. Indeed, although B is supported 
only over some small range of z centered around Zh/3, the metric functions A and S 
deviate from their vacuum values all the way from this point to the horizon, as required 
by causality. In contrast, the curvature outside the outgoing wave packet is very close to 
the curvature of the static black brane. 

These observations suggest that the states which take the longest time to thermalize 
are those that are initially localized close to the horizon on the initial-time slice. An 
example is provided by B{t = 0, z) ~ z^*^, whose evolution is shown in Fig. 2 (right). The 
reason is that the outgoing wave packet needs to escape the neighborhood of the horizon 
and travel all the way to the boundary to bounce off and finally fall into the black brane 
horizon. By localizing the initial profiles close to the horizon, the longest isotropization 
times that we were able to obtain with our numerics, which used rather moderate grids, 
were about 1.1/T — 1.2/T, with T the final equilibrium temperature (see Fig. 5). 

4 Linear approximation for holographic isotropization 

4.1 Leading order correction to the pressure anisotropy 

Linearizing Einstein's equations in the setup of holographic isotropization can be formally 
phrased as an expansion in the amplitude of perturbations on top of the AdS-Schwarzschild 
black brane. We thus write 

A(t, z) = \(l-z^) + a 5A^^\t, z) + Oia^), 

z'^ 

^{t,z) = - + a6Y.^'^\t,z) + 0{a^) and B{t,z) = a5B^^\t,z) + 0{a^), (4.1) 

where a is a formal parameter counting the order in the amplitude expansion. 

The smallness of the initial data can be physically quantified by either measuring the 
total entropy production on the event horizon or by following the amplitude of the pressure 
anisotropy during the evolution process and comparing it to the energy density. It is 
important to restress that we want to use the linearized approximation without necessarily 
restricting to the initial data being small perturbations of the AdS-Schwarzschild black 
brane, precisely in the spirit of the original close-limit approximation [23, 24] but now in 
the context of AdS spacetime. 

The initial data for the full nonlinear Einstein's equations are given by specifying the 
energy density £ and the form of -B as a function of the radial coordinate on the initial-time 
slice. As anticipated earlier, the main motivation for choosing B over S in specifying the 
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initial data was that the former appears quadratically in the constraint (3.1e). This feature 
persists also with the other components of the Einstein's equations apart from the equation 
(3.1b), which immediately leads to 

6A^^^t,z) = and 5sW(t,z) = 0. (4.2) 

6B^^'{t,z) on the other hand remains nontrivial and is a solution of the equation (3.1b) 
with A and S set to their form in the AdS-Schwarzschild background given in (2.9). 

The initial condition for this equation is the same as the initial condition for the full 
Einstein's equations, i.e. 

5B'^^\t = 0,z) = B{t = 0,z). (4.3) 

The energy density £, which is constant in our setup and is the remaining part of the initial 
state specification, is already included in the background we linearize on top of. 
In full detail, the equation for dB^^'{t, z) reads (with the choice of units £^ = |) 

;^(3 + z^) dJB^ -hi- z^) dl6B^^^ - ^ dt6B^^^ + d,dt SB^^^ = . (4.4) 

This is a simple, first order in time, partial differential equation which can be solved once 
supplemented with the Dirichlet boundary condition at z = 0. To improve the stability in 
our computations we actually worked with 

6Blll{t,z) = j,5B^'\t,z) (4.5) 

and required that 6Breg{t,z = 0) = 0. The other boundary of the computational grid is 
put inside the event horizon of the background solution. 

Equation (4.4) captures the leading order dynamics of the whole stress tensor, as 
the stress tensor is completely specified in terms of the energy density £ (encoded in the 
background) and the pressure anisotropy AV{t) (encoded in 6B^^'). 

In our previous work [1] we compared the predictions for AV{t) following from the 
nonlinear Einstein's equations with the ones obtained using (4.4) and we found that the 
results agreed with a 20% accuracy for the vast majority of the profiles we considered. The 
match of the qualitative features can be seen in Fig. 4. The histogram in Fig. 5 shows the 
quantitative match for the isotropization time for more than 800 non-equilibrium initial 
states which are all large perturbations of the AdS-Schwarzschild black brane. We see that 
the linearized Einstein's equations predict the isotropization time with a 20% accuracy in 
the vast majority of cases. 

An interesting feature visible in Fig. 4, as well as earlier in Fig. 3, is that the pressure 
anisotropy obtained from the linearized Einstein's equations always follows closely the full 
result at early times and, if at all, the curves start to differ only at some transient time. 
We understand this feature as a natural consequence of the fact that, due to the fixed 
asymptotics, the dynamics is approximately linear near the boundary of AdS, which is 
the bulk region causally responsible for the early-times dynamics. The pressure anisotropy 
curves start to differ only after the boundary is reached by a signal propagating from the 
interior of the geometry, as seen in Fig. 3, and this is the moment at which in the full 
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Figure 4. Comparison between the time evolution of the pressure anisotropy predicted by the linear 
equation (4.4) (red) and the full result (blue) for 15 different initial conditions. The leading order 
linearized Einstein's equations predict both qualitative and quantitative features of the dynamics 
of the dual stress tensor in our setup. A more thorough scan of the initial conditions (as shown in 
Fig. 5) did not reveal any instances in which the linearized approximation broke down. 

analysis nonlinear effects become most visible. Our main result can be phrased as the 
surprising statement that these linearities did not result in a large effect on the boundary 
stress tensor for any of the initial states that we analyzed. 
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Figure 5. (Top) Ai^^o is the difference between the isotropization time predicted by the full and 
the linear equations. The height of each bar in the histogram indicates the number of initial states 
for which the evolution yielded values in the corresponding bin. The total number of initial states 
is more than 800. We see both that holographic isotropization proceeds quickly, at most over a 
time scale set by the inverse temperature, and that the linearized Einstein's equations correctly 
reproduce the isotropization time with a 20% accuracy in most cases. Note that the histogram is 
based on a different sample of initial states than those originally considered in [1]. In particular, 
we incorporated the binary search algorithm absent in [1] and were stricter about the maximum 
violation of the constraint that we allowed. 

(Botom) Close inspection of one of the few profiles for which the linearized approximation seemingly 
fails by more than 20% (Atigo/tiso — —0.5) shows that it is the imperfect isotropization criterium 
which leads to the mismatch rather than the failure of the linear approximation. Indeed, the left 
plot shows that, on the scale of the initial anisotropy, the linear result yields a good approximation. 
However, the isotropization criterium makes no reference to this scale, and results in a 50% difference 
in the isotropization times, indicated by the arrows on the right plot. See [9] for a related discussion 
of subtleties involved in defining the thermalization (or more accurately hydrodynamization) time 
in a similar setup. 



16 



4.2 Connection with quasinormal modes 

Equation (4.4) can be solved either as an evolution equation given some initial profile for 
6B^^', as discussed in the previous section, or by decomposing 5B^^' as a superposition of 
modes with factorized time dependence: 



55W(t,z)~e-^'^^*6j-(z). 



(4.6) 



These modes are known as quasinormal modes, and they are characterized by the require- 
ments that they be normalizable near the boundary (z = 0) and that they obey ingoing 
boundary conditions at the event horizon {z = 1).^ The latter condition makes the frequen- 
cies coj complex with imaginary parts responsible for the exponential decay in time. The 
quasinormal modes (4.6) appear in pairs, as taking the complex conjugate of the equation 
(4.4) for the quasinormal mode with frequency ujj leads to the equation for the quasinormal 
mode with frequency — w*. This feature can be seen in Fig. 6. 
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Figure 6. The plot on the left shows the frequencies of the ten lowest quasinormal modes including 
their complex conjugates. The mode with the smallest negative imaginary part will by the dominant 
mode at late times. Notice that the spacing between the modes is approximately constant (it differs 
by about 0.1%). The plot on the right displays the lowest three quasinormal modes as a function of 
the radial coordinate z, where blue and red denote their real and imaginary parts. The normalization 
we use is such that the real part at the horizon {z/ Zh = 1) is equal to unity, whereas the imaginary 
part vanishes there. One clearly sees that higher modes (which decay faster) are more dominant 
near the boundary. 

In the context of gravitational collapse, the lowest quasinormal modes are known to 
govern the late-time decay of black hole perturbations (see e.g. [33]) and this is also ex- 
pected in the current setup. On the other hand, the results from [1], reviewed in the 
previous section, suggest that the equation (4.4) predicts rather well the full time depen- 
dence of the large-z behavior of the warp factor B. Hence a natural question is what is the 
quasinormal mode content of the perturbations that we considered. 

In order to answer this question we followed the prescription of [22] and computed the 
lowest 10 quasinormal modes (4.6) by solving equation (4.4) for the ansatz (4.6) in the 
near-horizon expansion and evaluating the resulting expression at the boundary to find 



^In the ingoing Eddington-Finkelstein coordinates the ingoing condition at the horizon is equivalent to 
regularity of the solution at the horizon [32] . 
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Uj's leading to normalizable modes. The (somewhat arbitrary) normahzation of our modes 
is fixed by demanding that at the horizon (z = 1) 

6,(1) = 1. (4.7) 

On Fig. 6 we plot the obtained frequencies lvj of the lowest 10 quasinormal modes, as 
well as bulk profiles for the real and imaginary parts of bi(z), 62(2:) and 63(2;) normalized 
according to (4.7). 

The idea now is to use the quasinormal modes to decompose solutions of (4.4), i.e. to 
write a solution of (4.4) in the form 



^B^QNMit,z) = Re 



-W^QNM 



i=l 



(4.8) 



where we truncated the expansion at some A'^qnm, although formally we could set A'^qnm = 00. 
In our calculations we used A'^qnm = 10. 

One can view (4.8) as a further simplification as compared with solving numerically 
(4.4), which approximates very well (within 20%) the full Einstein's equations when it 
comes to predicting the form of the dual stress tensor. The reason for this extra simplifi- 
cation is that now the solution is specified by providing a few complex numbers^ (say 10 
complex coefficients Cj's) which due to the linearity of the problem can be fitted on the 
initial-time slice to B(t = 0,z). 

As a way of generating coefficients Cj's we minimized 

'^\B{t = 0,z)-5B^^L{t = 0,z)\ (4.9) 

by using the least squares method on a discrete sample of the radial position z. Naturally, 
one needs far more points than the number of quasinormal modes included in (4.8). 

The subtlety in using (4.9) lies in the choice of the mutiplicative factor under the 
integral, which we set to be 1/z^. We checked that both 1/z and 1/z^ do not work well, 
as the first one does not take sufficiently into account and the other overcounts the near- 
boundary behavior of B{t = 0, z). On the other hand, l/z^ seems to work equally well as 
l/z"^, but for definiteness we focused here on the latter. 

Fig. 7 displays the difference between B{t = 0,z) and 6BQji^{t = 0,z) as a function 
of the number of quasinormal modes in two representative examples. Clearly, if a good fit 
is possible, then the profile (4.8) will solve the linearized Einstein's equations nicely since 
each quasinormal mode solves them individually. 

In Fig. 8 we compare the linearized evolution obtained from a direct solution of (4.4) 
and from a solution based on a decomposition into quasinormal modes. One can see 
that the contribution from each individual quasinormal mode can be large, ^ but that the 

*One may construct exceptional initial profiles, which are for instance very close to the boundary, or 
very rapidly oscillating. Including more quasinormal modes (taking A^qnm in (4.8) somewhat bigger than 
10) would allow us to treat these cases more accurately. 

^Note that this depends on the normalization of the individual modes. 
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final sum approximates the linearized evolution very well. Finally, in Fig. 9 we plot three 
representative examples, where the profile with B(t = 0, z) having support mostly in the 
IR displays this interference phenomenon particularly nicely. 
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Figure 7. The plot on the left displays the maximum of the error when approximating B{z) = 
B{z)/z^ by the first A^qnm (complex) quasinormal modes, with B{t = 0, z) = —ia^z^. The plot in 
the middle shows the error for the same profile as a function of the bulk coordinate z while using 
the 10 lowest quasinormal modes. The right plot displays the error for B{t = 0, z) = z^^ and shows 
clearly that a profile which is dominated in the IR is much harder to fit by the quasinormal modes. 
This causes oscillations in the evolution, as can be seen in Fig. 9. 
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Figure 8. On the left one sees the pressure anisotropy l^V /£ as predicted by the linearized 
evolution, or indistinguishably by the sum of the lowest 10 quasinormal modes as a thick blue line. 
One can also see there the sum of the first 1 (blue), 2 (green), 3 (orange) or 4 (red) quasinormal 
modes. As becomes apparent, the late time dynamics is well approximated already by keeping only 
the lowest quasinormal mode, but if one uses more the fit starts matching earlier. Note that the 
coefficients are computed such that the sum of the 10 fits the dynamics best. 

On the right we plot the individual quasinormal modes with the same coloring. One clearly sees 
that each of them carries very large anisotropy, but that their interference matches the linearized 
solution. 



4.3 Full leading order correction to the bulk 

As reported in [1] and briefly reviewed in section 4.1, linearizing the Einstein's equations 
leads to a very significant simplification in our setup when it comes to predicting the 
dynamics of the dual stress tensor, which is the only excited local operator. A natural 
question arising is whether also the dynamics of nonlocal observables, such as Wilson loops 
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Figure 9. In this figure we illustrate anisotropies of the full (blue), linearized (green) and quasi- 
normal mode (red) evolution of three representative initial profiles, located in the IR, spread-out 
and located in the UV respectively. Clearly, the initial profile located in the IR takes some time 
before exciting the anisotropy at the boundary, which also explains the late thermalization. The 
UV profile can have a very large anisotropy, but isotropises very fast. 

These features are nicely described when looking at the quasinormal modes coefficients Cj 's (below, 
real (blue) and imaginary (red) part). For the IR profile each individual contribution is very large, 
but they interfere in such a way to give only moderate anisotropies. In this way it is also possible to 
reach isotropization as late as 6 times the lowest QNM e-folding time. We also see that one would 
need to compute more quasinormal modes to accurately fit this profile. 

or the entanglement entropy, can be determined in a satisfactory way by using the linearized 
Einstein's equations alone. As nonlocal observables probing sufficiently large domain in 
a dual field theory obtain direct contributions from the interior of bulk spacetime, this 
question translates to which extend the linearized gravity reproduces the full bulk. 

Already at the superficial level it is easy to see that in order to answer this question 
we need to extend the analysis from [1] and include quadratic corrections to linearized 
Einstein's equations. A simple way to understand it is to notice that at the linear level 
both A and S retain their background form and hence the equation (4.4) does not capture 
the phenomenon of the entropy production. The latter appears only when A or Ti are 
different from their background values, as follows from equations (2.12) and (2.14). 

The lack of entropy production at first order in the amplitude expansion does not 
matter for initial states for which the initial condition B{t = 0,z) is tiny. In this case the 
quadratic corrections will be even tinier and as negligible as the actual entropy production 
in this particular process. However, we already saw that the linearized Einstein's equations 
do a good job in predicting AV{t) also in the cases when the perturbation is not small by 
any means. In these situations though, they would do a miserable job in predicting the 
entropy production. 

This argument implies that actually the full leading order correction to the bulk space- 
time in the current setup comprises 6B^^\ which captures the dynamics of AV{t) in the 
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leading order, and 6A^'^' and 6T,^'^' . The latter quantities are directly responsible for the 
leading order entropy production in the close-limit approximation in our setup. 

From a more general perspective, we will treat the entropy production as an example 
of a quantity that is highly sensitive to the IR part of the bulk geometry and we expect 
that the corrections SA^"^' and JS^^' will be crucial for computing other observables of that 
sort, including two-point functions, entanglement entropy, Wilson loops, etc. 

By looking closer at the full set of nonlinear Einstein's equations (3.1) one can easily 
anticipate that 6T,^'^> can be obtained by solving the constraint (3.1e) expanded up to the 
quadratic order in the amplitude 

za2 5S(2)+2 5,5S(2) = -l(a,J5W)2. (4.10) 

This is a second order linear ordinary differential equation in z, which, having obtained 
6B^^'{t, z), can be solved using very basic numerical methods on each time slice. The 
boundary condition for ST,^"^' is that it vanishes at z = 0, as the boundary behavior is 
always incorporated in the background. 

Having both B^^'{t,z) and T,^^>{t,z) we can use equation (3.1c) expanded to the 
quadratic order in amplitude to solve for 6A^'^'{t, z) 

z^d^SA^'^ + 2zdM^^ - 6^(2) = -^(1 - z'XdJB^'^f + Zdt5B^^^d-M^^ 

_1^(1 _ z4)(<5s(2) + zdM^^) + l2dt5T.^^\ (4.11) 

Again, one obtains a second order, linear, ordinary differential equation in z, which can be 
integrated with the use of basic numerical techniques. The choice of integration constants 
is dictated by the near-boundary expansion (2.6), which at low orders is the same for the 
linearized and the full Einstein's equations. 

At the second order of the amplitude expansion the correction to B vanishes. The 
reason is that equation (3.1b) at this order is homogeneous and first order in time, and we 
already used the initial profile for B as the initial condition for 5B^^' so that B^"^' {t = 0, z) = 
0. Finally, it is worth noting that in order to improve the accuracy of our procedures, in 
numerical calculations we actually redefined dT,^"^' and 5A^'^' to be 

(5S(| = z-55S(2) and 5 Af^ = z-H A'^''\ (4.12) 

Although formally we refer to our expansion scheme as 'the amplitude expansion', in the 
spirit of the close limit approximation we are actually interested in evaluating the whole 
leading order contribution to the geometry {6B^^\ 5T,^'^' and 5A"^) for initial data of 
arbitrary amplitude and in comparing it with the full answer. 

4.4 Leading order correction to entropy production 

Having obtained the geometry at leading order, i.e. having solved equations (4.4), (4.10) 
and (4.11), we use (2.12) and (2.14) to calculate the position of the event horizon as well 
as its entropy. 
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Contrary to the nonlinear case, we are not guaranteed that the entropy wiU be a non- 
decreasing function of time. An easy way to understand it is by noting that the solution 
of linearized equations corresponds to the solution of full equations with some 'matter' 
stress tensor in the bulk. As this matter is not guaranteed to obey the energy conditions, 
the approximated entropy may in principle decrease. Moreover, the linearized formula for 
the area of the event horizon following from e.g. (2.14) is not guaranteed to lead to a 
nonnegative result. 

It turns out, however, that most often these two caveats do not apply, and even if they 
do, they do not change the results significantly when compared with the full answer, as can 
be seen in Figs. 10 and 11, which are the entropy analogues of the anisotropy Figs. 4 and 
5. The results depicted in Figs. 10 and 11 show that the leading order correction consisting 
of 6B''^', 6Tj^'^' and SA^"^' reproduces a bulk IR-sensitive observable such as the entropy 
production with a 20% accuracy for the vast majority of states. The numerical techniques 
needed for obtaining these leading order approximation are much simpler than needed to 
obtain the nonlinear results. In particular, the built-in Mathematica command NDSolve 
is sufficient. This is already a significant simplification, as the actual procedures used to 
solve the nonlinear problem required a fully-fledged numerical implementation. 

4.5 Including higher order corrections 

Having computed the full leading order correction, a natural question is whether including 
higher order terms in the amplitude expansion will get us closer to the full nonlinear result. 
A direct calculation shows that the first subleading correction to the pressure anisotropy 
comes from 6B^^' , whereas the first subleading correction to the entropy production comes 
from the combined effect of including 6A^'^' and ST,^"^'. 

In order to understand the role of these subleading corrections, we compared the 
predictions for the pressure anisotropy from the following three approaches: 

• nonlinear simulations; 

• leading order solution of linearized equations, i.e. 6B^^\ (5S'^' and 6A^'^'; 

• solution of the linearized equations including the first subleading corrections, i.e. 6B''^>, 
(5S(4) and 5^(4) . 

We started with some initial B{t = 0, z) and, after each run, we increased its amplitude by 

fixed amount. For definiteness, in Figs. 12 and 13 we show the results for profiles obtained 

from 

B(t = 0,z) = -£z^ (4.13) 

6 

by increasing the amplitude by ^£ at each run. As expected, the leading order approxi- 
mation predicts well both the pressure anisotropy and the entropy production for all am- 
plitudes, including the largest ones that we could simulate. Also unsurprisingly, for small 
enough amplitudes the subleading correction improves the result as compared to the lead- 
ing order result. However, for sufficiently large amplitudes the expansion does not seem to 
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Figure 10. Comparison between the leading order prediction for the entropy production coming 
from the Einstein's equations expanded to quadratic order (red curves) and the exact entropy 
production the obtained from the fuU nonhnear analysis (blue curves), for the same 15 initial states 
as those in Fig. 4 with the same ordering. Note that none of the initial conditions can be regarded 
as a small perturbation either in the sense of the amplitude or the total entropy produced during 
the equilibration process. On Fig. 11 it can be seen that the leading order prediction for the entropy 
production from linearized Einstein's equations matches remarkably well, within 20%, with the full 
nonlinear result for the vast majority of the initial states that we analyzed. 

converge. In the case of the pressure anisotropy, including the subleading correction over- 
shoots the prediction for the initial states with the largest amplitudes. In the case of the 
entropy density, including the subleading correction for large-amplitude profiles leads to 
rather unphysical results since the entropy density decreases over certain periods of time. 
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Figure 11. Histogram for the error in the entropy production analogous to that in Fig. 5 for 
the isotropization time. Both histograms are based on the same sample of more than 800 ini- 
tial states. We define the error in the entropy production as AEntropy — a7 ' where 
As — s{t — cxd) — s{t = 0) and As^^^ = s(^)(i = cxd) — s^^\t = 0) are the exact and the leading 
order entropy productions, respectively. We see that the leading order approximation reproduces 
the correct result with a 20% accuracy for the vast majority of states. 

A further indication that adding subleading corrections does not improve the result for 
generic states conies from Fig. 14, where we see that the histogram for Atis^/tiso does not 
become narrower along the Ati^o axis after including the third order correction for A'P(t). 
With the prospect of studying other more complicated setups, all these results point to 
the conclusion that it is safe, and generically better, to keep only the leading order terms. 

5 Summary and open problems 

In this article we presented a thorough study of the simplest example of holographic ther- 
malization process, the holographic isotropization. The primary motivation for these inves- 
tigations was understanding to which extent the full process of holographic isotropization 
can be simplified by describing it using the linearized Einstein's equations alone. Such an 
approximation is an AdS generalization of the closed limit approximation for asymptoti- 
cally flat spacetimes, in which in certain black hole mergers one linearizes the Einstein's 
equations on top of the final black hole predicting quite accurately the gravitational wave 
pattern at infinity. Our previous studies in [1] showed that such generalization of the close 
limit approximation works remarkably well, i.e. with a 20% accuracy, when it comes to 
predicting the time dependence of the expectation value of the boundary stress tensor. 
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Figure 12. Pressure anisotropy as a function of time as predicted by the full Einstein's equations 
(blue curves), by the first order approximation (red curves) and by the first +third order approx- 
imation (green curves). In all cases the initial state is of the form (4.13) with an amplitude that 
increases by |£ from one plot to the next (from left to right and from top to bottom). We see that 
for small (large) amplitudes adding the third order result improves (worsens) the agreement with 
the exact result. 
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Figure 13. Entropy of the event horizon as a function of time as predicted by the full Einstein's 
equations (blue curves) , by the second order approximation (red curves) and by the second+fourth 
order approximation (green curves). In all cases the initial state is of the form (4.13) with an 
amplitude that increases by :^8 from one plot to the next (from left to right and from top to 
bottom). We see that the second order prediction works better for all cases shown. 

The main phenomenological motivation for studying holographic thermalization is 
learning possible lessons about the way the thermalization (or rather hydrodynamization) 
process proceeds in relativistic heavy ion collisions at RHIC and LHC. By replacing QCD 
by a theory with a gravity dual one only expects to obtain either qualitative insights or 
quantitative ball-park estimates [34]. In this sense a 20% accuracy is more than what is 
needed in order to understand the phenomena we are interested in, and at the same time 
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Figure 14. Same histogram as in Fig. 5 except that the first order approximation for the pressure 
anisotropy is replaced by the first+third order approximation. Comparing with Fig. 5 we see 
that the histogram does not become narrower along the Ati^o axis after including the third order 
correction. 

may allow to address otherwise technically hard-to-tackle questions. Two examples of such 
problems in the relativistic heavy ion collisions context are the pre-equilibrium phase of 
the elliptic flow and the equilibration of transverse-plane inhomogeneities following from 
event-by-event fluctuations. Solving their holographic analogues in full generality will re- 
quire complicated simulations of AdS-black hole spacetimes depending on all five bulk 
coordinates and our hope is that a suitably developed linear approximation may allow us 
to obtain results with a reasonable accuracy at a much smaller cost. 

In addition to elaborating in detail on our previous studies, the chief aim of the current 
paper was to establish to what extent our approximation reproduces the full thermalization 
process, i.e. not only the behavior of fields in the vicinity of the boundary (responsible for 
the expectation values of gauge theory operators), but also the entire bulk spacetime (which 
carries direct information about e.g. nonlocal probes in the gauge theory). 

In our previous studies we linearized Einstein's equations around the final-state black 
brane solution and we kept only first order terms in the amplitude expansion. We then 
solved these linear equations supplemented with appropriate AdS boundary conditions and 
with the same initial conditions as required to solve the full Einstein's equations. In the 
spirit of the closed-limit approximation the initial conditions were not necessarily small 
perturbations of the final black brane. This procedure correctly captured the evolution 
of the expectation value of the gauge theory stress tensor with a 20% accuracy, but it 
predicted no entropy production despite the fact that this was often substantial. 

From this perspective, our previous studies in [1] provided only a partial contribution 
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to the full leading order approximation to the bulk spacetime, which also includes terms 
quadratic in the amplitude expansion. These quadratic terms do not contribute to the 
expectation value of the stress tensor, but they do modify the bulk geometry. In order to 
probe this effect, in this paper we considered the contribution of the quadratic terms to 
the entropy production and concluded that this is also reproduced with a 20% accuracy. 

We also analyzed the effect of higher order corrections in the amplitude expansion. Not 
surprisingly, we found that for initial states that start off sufficiently close to equilibrium the 
inclusion of these terms improves the agreement with the nonlinear evolution. However, for 
large deviations from equilibrium (large amplitudes) the inclusion of these terms actually 
spoils the 20% agreement between the leading order approximation and the exact result. 
This suggests that in future cases in which a full nonlinear simulation might not be feasible, 
it might be enough to work with the leading order approximation. 

An interesting byproduct of our analysis is that the quasinormal modes, which so 
far were thought to be responsible only for the late time approach to equilibrium in the 
expectation values of dual operators, might actually predict quite accurately the full time 
dependence of the dual stress tensor provided a sufficient number of them is included. 

The most important open problem raised by our previous work [1] and the current 
paper is how to generalize the close-limit approximation to more phenomenologically inter- 
esting setups, for example to expanding plasmas such as those considered in [6-8]. Those 
setups are qualitatively different form the one we have considered because of two interre- 
lated features. First, the late time physics will be described by a solution of hydrodynamic 
equations which is not known a priori in terms of the initial conditions without solving for 
the evolution of the system. This contrasts with the current case, in which the final state 
is a homogeneous thermal state whose energy density is known from the start. Second, 
the late time state will not be static, so one does not expect to be able to describe it by 
linearizing around the AdS-Schwarzschild black brane. Instead, since the energy density 
of an expanding plasma gets more and more diluted as time progresses, one might naively 
hope to be able to linearize around the Poincare patch of vacuum AdS. An indication that 
this is not the right procedure comes from Ref. [21], whose authors analyzed (mainly) the 
gravitational collapse of a massless scalar field triggered by turning on its non-normalizable 
mode with a small amplitude. The naive small amplitude expansion on top of vacuum AdS 
had to be resummed so that the background became the AdS-Vaidya solution instead. The 
latter geometry has a macroscopic horizon which governs the dissipation of perturbations 
propagating on top of it. This suggests that, in cases of expanding plasmas, the correct 
procedure is to linearize around some other expanding configuration. We will report else- 
where on studies in this direction in the context of the simplest expanding plasma system, 
the boost-invariant flow [35]. 
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